function d = loadDisturbance(dt,n,rand)
% This does the same thing as makeLoadProfile but then interpolates d to
% sample at a rate of dt for n zones. Set rand = 1 for random load, and set
% rand = 0 for checkign consistency (will use first n profiles).
% makeLoadProfile:
% Take an average of load data over one day from Doe Library data and put it in 
% Pd.time and Pd.values as well as Toa.time and Toa.values.
% Each row in Pd indicates a different zone's averaged load.

%% MPC Lab Data
% get MPCLab.load, MPCLab.time where each row is a different day
load MPCLabLoad

t0 = MPCLab.time(1,1);
tend = MPCLab.time(end,end);
numdays = tend-t0;
t = t0:dt:(t0+1);

% randDays = ceil(rand(1,n)*(numdays-1));
orderedDays = 1:n; randDays = orderedDays;

d = zeros(n,length(t));

% map a zone to a random day in MPC lab data
for zone = 1:n
    day = randDays(zone);
    t0zone = MPCLab.time(day,1);
    d(zone,:) = interp1(MPCLab.time(day,:),MPCLab.load(day,:),t-t0+t0zone,'linear','extrap'); 
end

%% Debug MPC lab by plotting it
% close all
% figure(10)
% for zone = 1:n
%     plot(t,d(zone,:))
%     hold on
% end
% title(['Disturbance for ' num2str(n) ' zones'])
% datetick('x','HH');
% ylabel('d (^o C)')
% xlabel('time (HH)');

%% HW 7 Data
% % Using HW7_data, this gave us 5 minute sampled data on July 3, 2009
% % load Pd, T_oa
% load HW7_data
% Pd.time = repmat(Pd.time,zones(end),1);
% Pd.values = repmat(Pd.values',zones(end),1);
% T_oa.values = T_oa.values';
% T_oa.values = T_oa.values-273.15; % Convert to degC

%% DOE Library Data
% load DOE_data.mat
% % Format of DOE library data:
% % Each cell is a different zone.
% % Each cell contains different data structs, each struct contains a matrix
% % time and a matrix data. The ROW corresponds to DAY and the COLUMN
% % corresponds to time in that same day.
% 
% % Let's take data from the first and second zones and average it.
% zones = 1:8;
% % deltat = 15*60;
% deltat = 1;
% [numDays,numTimePoints] = size(DOE_data(1,1).LoadEst.time);
% Pd.values = zeros(zones(end),numTimePoints+1);
% T_oa.values = zeros(1,numTimePoints+1);
% Pd.time = [DOE_data(1,1).LoadEst.time(1,:),DOE_data(1,1).LoadEst.time(1+1,1)];
% T_oa.time = Pd.time;
% Pd.time = repmat(Pd.time,zones(end),1);
% for i = 1:numDays-1
%     for zone = zones
%         p = DOE_data(1,zone).parameter;
%         % We have to append the first data point on the next day to get a full
%         % 24 hours worth of data
%         Pdtemp_d = [deltat*DOE_data(1,zone).LoadEst.data(i,:)+p(1)*DOE_data(1,zone).Toa.data(i,:)+p(5)*DOE_data(1,zone).SunRadiation.data(i,:),...
%                     deltat*DOE_data(1,zone).LoadEst.data(i+1,1)+p(1)*DOE_data(1,zone).Toa.data(i+1,1)+p(5)*DOE_data(1,zone).SunRadiation.data(i+1,1)];
%         % Divide by (numDays-1) = number of averaging points and add
%         Pd.values(zone,:) = Pd.values(zone,:) + Pdtemp_d./(numDays-1);
%     end
%     
%     % Collect Outside Air Temperature
%     Toatemp_d = [DOE_data(1,zone).Toa.data(i,:),DOE_data(1,zone).Toa.data(i+1,1)];
%     T_oa.values = T_oa.values + Toatemp_d./(numDays-1);
% end
% 
% % Convert T_oa from Kelvin to Celcius
% T_oa.values = T_oa.values - 273.15;

%% Debug
% close all
% cc = hsv(12);
% % zones = 1:8; %IF using DOE library
% for zone = zones
%     plot(Pd.time(zone,:),Pd.values(zone,:),'color',cc(zone,:));
%     hold on
% end
% legend('1','2','3','4','5','6','7','8')
% title('Load data');
% figure
% plot(T_oa.time(1,:),T_oa.values(1,:));
% title('Outside Air Temperature Data');

%% Set up discretization at correct sampling rate
% t0 = Pd.time(1,1);
% tend = Pd.time(1,end);
% 
% C = 9.2e3*eye(n); % kJ/degC 
% R = 50*eye(n); % degC/kW 
% 
% % determine disturbance for MPC controller sampling
% datevec = t0:dt:tend; 
% d = zeros(n,length(datevec));
% for i = 1:length(datevec)
%     t = datevec(i);
%     d(:,i) = dt*C^-1*ones(n,1)*interp1(Pd.time(3,:),Pd.values(3,:),t)+...
%             dt*R^-1*C^-1*ones(n,1)*interp1(T_oa.time,T_oa.values,t);
%     d(1,i) = dt*C(1,1)^-1*interp1(Pd.time(4,:),Pd.values(4,:),t)+...
%             dt*R(1,1)^-1*C(1,1)^-1*interp1(T_oa.time,T_oa.values,t);
% end